import os
import matplotlib.pyplot as plt
import numpy as np
import numpy.matlib
import math

os.chdir('D:\\ProjetQ\\methode_estimation')

global Ptplan #for recording all the values of the chosen plan 
global Ptvaleur
global echelleX,echelleY 
global npas, pas, ori

ori=range(3)
npas=range(3) #npas : numbers in each direction ,y,z
pas=range(3) #pas : distance in each direction x,y,z

f1="UreaDens.cube"
f2="pt_plan.txt"

##################################################
#1 general info  
def readtitle(fichier): #title
    fichier.readline()
    fichier.readline()  #test function

def infoAtomes(fichier): #atom information reading
    v=fichier.readline()
    v=v.split()
    Nt=int(v[0]) #total atom numbers
    ## print 'The total atom numbers in this molecule is', Nt

    for i in range(3): #origin information reading
        ori[i]=float(v[i+1])
    ## print 'The l\'origin postition of the molecule is', ori

    #step number reading
    for i in range(3):
        v=fichier.readline()
        v=v.split()
        npas[i]=int(v[0])
        pas[i]=float(v[i+1])
    
    for i in range(Nt): #jump over the atom information for the moment
        fichier.readline()


#################################################
#2 value of a single point in the .cube
        
def bloc(Nz): #line numbers for a block, Nz=npas[2]
    t=Nz/6
    if Nz%6:
        t=t+1
    return t

def sauterbloc(fichier,nx,ny):
    Nbloc=nx*npas[1]+ny #For the point(nx,ny,nz), the number of blocks to be jumped over
    LnInBloc= bloc(npas[2]) #line number in one block
    for i in range(Nbloc):
        for j in range(LnInBloc):
            fichier.readline()
    
        
def ptori(fichier,nx,ny,nz):#find the value of a specific position
    fI=open(fichier,"r")
    readtitle(fI)
    infoAtomes(fI)
    
    if not -1<nx<npas[0]:
        raise ValueError, 'The defined area of x is depassed.'
    if not -1<ny<npas[1]:
        raise ValueError, 'The defined area of y is depassed.'
    if not -1<nz<npas[2]:
        raise ValueError, 'The defined area of z is depassed.'

    sauterbloc(fI,nx,ny)

    nz=nz+1
    Vertical=nz/6 #read the bloc which contians the needed information 
    Horizon=nz%6
    if not Horizon:
        for i in range(Vertical-1):
            fI.readline()
        Horizon=6
    else:
        for i in range(Vertical):
            fI.readline()
    v=fI.readline()
    v=v.split()
    valeur=float(v[Horizon-1])
    return valeur
    fI.close()

def ptsearch(x1,x2,x3): #search for the value that is not mentioned in the file
    point=[x1,x2,x3]
    rang=range(3)
    for i in range(3):
        for j in range(npas[i]-1):
            if j<=point[i] and (j+1)>point[i]:
                rang[i]=j
                break
 
    val1=ptori(f1,rang[0],rang[1],rang[2])
    val2=ptori(f1,rang[0],rang[1],rang[2]+1)

    p1=x3-rang[2]
    p2=1-p1
    valeurFin=p1*val1+p2*val2
    print 'valeurFin',valeurFin
    return valeurFin
            
    
##########################################
#3 (discretement)locate the section passing through 3 determined points in the grid, search the position on the axe z

def dotproduct(X1,X2):
    x=X1[1]*X2[2]-X1[2]*X2[1]
    y=X1[2]*X2[0]-X1[0]*X2[2]
    z=X1[0]*X2[1]-X1[1]*X2[0]
    return np.array([x,y,z])

def zerotest(X): # to test if a vector is zero
    d=0
    for s in X:
        d=d+abs(s)
    return d

def plan(fichier,X1,X2,X3): #X1,X2,X3 are vectors in 3D, this function is to search for all the positions and their values of the plan
                            #that we have defined. We only consider the case of the axis Z in this function
    Ptplan=np.matlib.zeros((npas[0],npas[1]),dtype=float)
    Ptvaleur=np.matlib.zeros((npas[0],npas[1]),dtype=float)    
    
    r=[X1,X2,X3]
    oriplan=np.array(X1)-np.array(ori) #We choose X1 as the origin of the plan
    V1=np.array(X2)-np.array(ori)
    V2=np.array(X3)-np.array(ori)

    Npas=np.array(pas) #normalize the positions in the formalisme of part 2 (?)
    oriplan=oriplan/Npas
    V1=V1/Npas
    V2=V2/Npas
    
    Vplan=dotproduct(V1,V2) #Vplan*Ptplan=Vplan*oriplan=ValeurPlan

    ValeurPlan=0
    temp=Vplan*oriplan
    for i in temp:
        ValeurPlan+=i

    for i in range(npas[0]):
        for j in range(npas[1]):
            Ptplan[i,j]=(ValeurPlan-i*Vplan[0]-j*Vplan[1])/Vplan[2]

    echelleX=pas[0]*math.sqrt(1+(Ptplan[1,0]-Ptplan[0,0])**2) #real
    echelleY=pas[1]*math.sqrt(1+(Ptplan[0,1]-Ptplan[0,0])**2)

    fO=open(fichier,'w')
    string=str([echelleX,echelleY])
    fO.write(string)
    string=str([npas[0],npas[1]])
    fO.write(string)
    
    for i in range(npas[0]):
        for j in range(npas[1]):
            Ptvaleur[i,j]=ptsearch(i,j,Ptplan[i,j])
            string=str([i,j,Ptvaleur[i,j]])
            print 'str',string
            fO.write(string)
    fO.close()
        
    return
##################################################
#Part trace
def sectrace(matr,dspx,dspy,nspx,nspy,xori,yori):
    plt.figure(1)    # create the figure

    # search the min and max of values in the matrix
    maxi=mini=matr[0,0]  #initialize maxi and mini
    for j in matr:
        if j>maxi:
            maxi=j
        if j<mini:
            mini=j
    V=np.arange(mini,maxi,(maxi-mini)/100)  # define the value of each contour
    X=[xori,xori+(nspx+1)*dspx,dspx] # coordinate the x axis
    Y=[yori,yori+(nspy+1)*dspy,dspy] # coordinate the y axis
    CS=plt.contour(X,Y,matr,V,linewidths=2)
    plt.show()
    return

##################################################
#Part test

ftitle=open(f1,'r')
readtitle(ftitle)
infoAtomes(ftitle)
ftitle.close()

X1=[0,0,1]
X2=[0,0,0]
X3=[0,1,0]
plan(f2,X1,X2,X3)

sectrace(Ptvaleur,echelleX,echelleY,npas[0],npas[1],0,0)

